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ABSTRACT 

We calculate the structure of a force-free magnetosphere which is assumed to coro- 
tate with a central star and which interacts with an embedded differentially rotating 
accretion disc. The magnetic and rotation axes are aligned and the stellar field is 
assumed to be a dipole. We concentrate on the case when the amount of field line 
twisting through the disc-magnetosphere interaction is large and consider different 
outer boun dary conditions. In general the field line twisting produces field line in- 
flation (eg. Bardou fc Heyvaerts 1996 ) and in some cases with large twisting many 
field lines can become open. We calculate the spin-down torque acting between the 
star and the disc and we find that it decreases significantly for cases with large field 
line twisting. This suggests that the oscillating torques observed for some accreting 
neutron stars could be due to the magnetosphere varying between states with low and 
high field line inflation. Calculations of the spin evolution of T Tauri stars may also 
have to be revised in light of the significant effect that field line twisting has on the 
magnetic torque resulting from star-disc interactions. 

Key words: accretion, accretion discs - stars: magnetic fields - stars: rotation - 
MHD - stars: neutron - stars: pre-main-sequence. 



1 INTRODUCTION 

Accreting magnetised stars exist in the form of neutron stars in close binary systems and T Tauri stars with surrounding 
protostellar discs. In both cases an accretion disc is disrupted through interaction with the stellar magnetosphere and the 
accretion flow becomes channelled along field lines. Torques act between the star and the accreting material and they determine 
the spin history of the star and its equilibrium rotation period. A spin-up torque is produced by rapidly rotating material in 
the inner regions of the disc while a spin-down torque is produced by the more slowly rotating material in the outer parts. 
Calculation of the torques requires knowledge of the structure of the disc and the magnetosphere in which it is embedded. 

Early work ( |Ghosh fc Lamb 1979a : Ghosh & Lamb 1979b ) suggested a direct dependence of the total torque on the mass 
accretion rate. However, recent observations of accreting neutron stars (Nelson et al. 1997; Shakrabatry et al. 1997) suggest 
that the torques can oscillate producing alternat ing spin-up and spin-down p hases with no evidence of a consistent correlation 



between the torque and the mass accretion rate. Li fc Wickramasinghe (1998 ) suggest that this may be due to a variable outer 



magnetosphere which causes a corresponding variation in the spin-down torque acting between the star and the disc. In this 
paper we investigate the structure of a magnetosphere corotating with the central star. We calculate the spin-down torque 
due to the star-disc interaction and we find that it is very sensitive to the amount of field line twisting that occurs through 
the coupling of the stellar field lines to the disc. The magnetospheric structure and the resulting torque are also sensitive to 
the applied boundary conditions. In the present work we neglect any torque that may arise from a stellar wind. 

According to the Ghosh & Lamb model the magnetic field of the central star penetrates the accretion disc due to 
diffusion arising from the growth of various instabilities. In the outer parts of the disc, viscous stresses are more important 
than magnetic stresses so that the effects of the magnetic field on the disc structure can be ignored. In the inner parts, the 
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flow becomes dominated by magnetic stresses which eventually lead to the truncation of the disc at an inner radius, 7?d, and 
the channeling of the flow to the star along stellar magnetic field lines. The twisting of the poloidal field by the vertical shear 
is assumed to be counteracted by the reconnection of the oppositely directed toroidal fields through the disc vertical extent. 
This occurs on a rapid timescale resulting to an equilibrium ratio for toroidal, B v , over vertical, B z field components at the 
disc surface. Alternatively the growth of B v can be controlled by turbulent diffusion through the disc ( Campbell 198?| ; Jy] 
1994|)]or reconnection of field lines in the magnetosphere (i.e. Livio & Pringle 1992). In all these cases magnetic torques of a 
similar form are produced (Wang 1995). 

In all calculations of the magnetic torques resulting from disc-star interactions B z is often assumed to be approximately 
equal to the stellar dipole field. However, twisting magnetic field lines imply the existence of currents which will be distributed 
in the magnetically dominated magnetosphere which is assumed to corotate with the central star. The star-disc magnetosphere 
is expected to attain a force-free equilibrium in an Alfven wave crossing time, as is the case in the solar corona. In the force- 
free equilibrium, the poloidal magnetic field will deviate from that of a dipole, the effect being especially significant in the 
case of large twisting (or \B V /B Z \). We find that the magnetic torque acting on the disc may also be significantly reduced. 

The related problem of the evolution of an initially potential field in response to an increasing footpoint shear through 
a sequence of quasi-static force-free equilibria has been addressed in the solar context. But note that, unlike in the work 
presented here, ideal MHD is assumed. Numerical calculations supported theoretical expectations (Aly 1984; Aly 1985; Aly 
1991 ; [Sturrock 1991 ) that the energy of the force-free field increases monotonically with increasing shear, approaching the 
energy of an open field (Klimchuck 1991; Roumeliotis, Sturrock & Antiochos 1994; Mikic & Linker 1994). Analytical sequences 



of quasi-static equilibria were also constructed using a self-similar approach (Newman, Newman & Lovelace 1992; Lynden 



Bell & ; Boily 1994; Wolfson 1995) and it was shown that the complete opening of the initially closed field lines is possible at 
a finite shear. This is associated with the development of current sheets as was conjectured by Aly (1985). 

Bardou & Heyvaerts (1996) have calculated magnetospheric force-free equilibria resulting from the twisting of an initially 
dipolar field anchored on a star which threads an embedded differentially rotating disc. Disc resistivity was taken to be of 
turbulent origin. However, because of problems associated with their numerical method these authors were not able to calculate 
equilibrium configurations for \B V /B Z \ larger than 2.55 at the surface of the disc. In most cases a maximum \B V / B z \ ~ 1 was 
adopted. Nonetheless significant field line inflation was observed. 

In the work presented here, we similarly find that moderate to large shearing of the stellar dipole field through interaction 
with the disc can lead to a partially open configuration, with expulsion of toroidal field, for favourable outer boundary 
conditions. Although the process can be inhibited by the presence of an outer boundary that confines the field, significant 
field line inflation may occur with a large reduction of the magnitude of the vertical field in comparison to the initial dipole 
value. This could h ave important consequences for studies of stellar sp i n evolut i on which have assumed a dipolar pol oidal 
field at equilibrium ( |Ghosh fe Lamb 1979b| ; |Cameron fc Campbell 1993[ |Yi 1995| ; |Ghosh 1995| ; |Armitage fe Clarke 1996; ). 

After formulating the physical model and basic equations in Section |2j, we describe our numerical method in Section ^| 
Using this we were able to perform calculations with a maximum value of \B V /B Z \ at the disc surface of 120 and with a 
variety of external disc radii. Results for Keplerian discs with different amounts of field twisting are given in sections ^. These 
calculations were performed with a field-confining, perfectly conducting outer boundary. 

We also per formed calculations with a partly sub-Keplerian accretion disc using the form of angular velocity, Q given by 
Campbell (1987 ). T his emulates the existence of an inner boundary layer as proposed in the Ghosh & Lamb model. Results are 
given in Section L3 for the conducting outer boundary condition. In Section we investigate the effect of the outer boundary 
condition, showing that a boundary condition that allows field lines to penetrate it, leads to more open configurations. 

In Section 4.5 we go on to evaluate the spin down torque acting between star and disc for the force free configurations. 
We show that this torque can be reduced by up to a factor of 100, in comparison to what would be obtained assuming an 
unmodified dipole field, in a configuration that has undergone large twisting and toroidal field expulsion. 
Finally in Section |^ we summarize and discuss our results and their application. 



2 BASIC EQUATIONS 
2.1 Force— free equilibria 

We consider a magnetic field £?d P anchored on a star that rotates with a rate fi*. We suppose that the magnetic field is 
axisymmetric with symmetry axis aligned with the rotation axis. A thin accretion disc orbits in the equatorial plane such 
that the system remains axisymmetric. We ignore the internal dynamics of the disc and assume that the material rotates 
with a prescribed angular velocity £1. For the most part this will equal or be close to the local Keplerian angular velocity 
Qk = \J GM/R 3 where R is the radial distance to the star measured in the midplane of the disc. We suppose that there is a 
highly conducting low density corona that corotates with the star and in which the disc is immersed. 

However, the disc has non zero resistivity such that the stellar field is enabled to diffuse into the disc. Because the 
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star/corona and disc rotate at different rates, field lines that permeate the disc are twisted in the azimuthal direction such 
that a toroidal component of the magnetic field is generated. 

The disc is considered to be turbulent and manifesting an effective magnetic diffusivity, r\. This enables a steady state to 
be achieved in which the generation of toroidal magnetic field through differential rotation between the disc and star/corona 
is balanced by the effects of turbulent diffusion. 

The corona responds to the twisting of the initially dipolar field by producing toroidal magnetic field and current density 
components. These adjust such that a force-free equilibrium is eventually attained. We expect this situation to occur when 
the coronal Alfven speed sufficiently exceeds the characteristic rotational speeds. This requires a low density corona which 
is magnetically dominated. For the purposes of the work presented here, we assume that such a corona exists in which a 
force-free equilibrium is established in a time short compared to the stellar rotation period. Thus we look for steady state 
axisymmetric configurations involving the star, the disc and the corona. 

An axisymmetric magnetic field B can be split into poloidal B p and toroidal B v ip components such that: 

B p = Vx(*y/(rsin60), (1) 

where is the magnetic stream function which is constant on field lines. Here we use spherical polar coordinates (r, 9, ip) 
based on the central star. Then 

_ < 1 d ® i afr- 

p V r 2 sin 6 d6 ' r sin 9r . 

The force-free condition implies that the Lorentz force is zero everywhere, thus: 



(2) 



JxB = (3) 

where J = (c/47t)V x B is the current density (c.g.s units are used throughout). The toroidal component of equation Q 
gives: 

B p ■ V(rsin6»B ¥ ,) = 0. (4) 

The definition of ^ implies constancy on field lines or equivalently B p ■ V\& = 0. Thus (Q) implies that: 

rsm6B v = /(*) = /, (5) 

where / is an arbitrary function of $ . 

Using (^) and the poloidal component of (^) we finally obtain the governing equation for in the form: 

A* = -//' (6) 
where /' = d//d\l/ and A is the differential operator defined through: 

A- — + — i^— (7) 
dr 2 r 2 3/i 2 ' 

where /i = cos 6. 

Equation (|^) can be solved for ^ if / is known. Since f(^) is constant on poloidal field lines we need to calculate its value at 
only one point along each field line. This we do by consideration of the interaction of the field with the resistive disc. 



2.2 Disc— stellar field interaction 

The corona above the disc is assumed to corotate with the star and so a toroidal field will be generated due to the vertical 
shear. Neglecting the radial component of the field interior to the disc, we calculate the internal toroidal field using the toroidal 
component of the induction equation which in cylindrical coordinates (R, ip, z) takes the form: 

^ = RB P • Vfi + ^A C (RB V ) + -3-V77 • V(RB ¥ ) (8) 
ot rt it 

where A c = V 2 — (2/ R)(d/dR) and 77 = T>v is the turbulent magnetic diffusivity which is assumed to be proportional to the 
turbulent viscosity v, the constant of proportionality being T>. We use the Shakura & Sunyaev (1973) parameterization for u, 
v = ass^K-ff 2 . Here the viscosity parameter is ass and H is the disc semithickness. The gas velocity in the disc is taken to 
be u = (0,rfi,0). 

We consider that B v is generated by the shearing of the vertical magnetic field component arising from the z dependence 



of f2. For simplicity we take 77 to be a function of radial distance only. For a more general treatment see Campbell (1987). 
Noting that the disc is thin we expect d/dz ~ (R/ H)d /dR and so we neglect radial derivatives in equation (^). Then in a 
steady state we have: 

d 2 B v m 

ri-^r = -RB Z — (9) 
oz z oz 
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To find B v we need to specify dQ/dz. We follow previous authors (Wang 1987; Campbell 1987; Yi 1994) and take the vertical 
shear to be concentrated in a thin layer close to the disc surface. There the angular velocity changes between the disc midplane 
value and £7*. We take the vertical average of equation the right hand of which gives: 



RB Z 
H 



dz 



i?B z (ft* - fi) 
H : 



(10) 



where B z is assumed not to vary with z thus being the disc midplane value. In the following D. will refer to the disc midplane 
value. 

We assume a simple approximation for the vertical average of the dissipative term on the left-hand-side of equation (^) namely: 

rll 



H J V 3^ 



dz 



-7/7 



B v 

JJ2' 



(11) 



where B v now applies to the value at the upper disc surface and 7 is a dimensionless parameter expected to be of order unity. 

Combining equation ( E^ ) with equation ( J]~l| ) we obtain the following estimate for the equilibrium value of the ratio B v jB z 
at the upper disc surface: 



?1 



77/ 



(12) 



We have inserted the ± alternative to allow for different senses of rotation of the star and disc with respect to the right 
handed coordinate system with the negative sign corresponding to clockwise rotation. The angular velocities are then always 
positive. Using the turbulent diffusivity specified above and taking 7 = 1, equation ( |l2j ) gives for Q — Qk' 

R\ 3 / 2 



B z 



± 



R 

Va ss H 



(I)' 



(13) 



From equation (^|) we see that B v changes sign at the corotation radius, R c = (GM/Q 2 ) 1 ^ 3 . Similarly the vertically 
integrated zip component of the magnetic torque per unit area, B Z B V R 2 /(2n), also changes sign such that the star transfers 
angular momentum to the disc outside the corotation radius while it gains angular momentum from disc material inside that 
radius. 

In our calculations we shall assume that RHQ,*/ (77/) = C = constant. Equation (^) then becomes (adopting the negative 
sign corresponding to clockwise rotation): 

3/2 



T z - C 



(I)' 



(14) 



Our assumption of C = constant is such that for large values of R/R c , B v / B z « — C which is constant. However, an important 
result of the calculations, namely the inflation of poloidal field lines such that B z — * at the disc surface for large C, should 
not depend on having a dependence of C on radius provided it remains large. For H/R — 0.1, ass = 0.01 and T> = 7 — 1 at 



the corotation radius, these being reasonable values for accretion discs, equation ( |l4"|) gives \B V /B Z 
We comment that equation ( |l3| ) indicates that in a disc with constant H/R and constant ass, 



Coc 



V \Ro) 



3/2 



C = 10 s for i? » i? c 



(15) 



Thus in this case if T> = 1, C would indeed increase with radius. 

Our analysis is distinct from that of e.g. Livio & Pringle (1992) who use a similar expression for \B V IB Z \ to equation (^J) 
(for R ^ R c ) but who assume that fast reconnection in the star-disc corona with the coronal Alfven speed restricts \B v jB z \ 
to a maximum value of unity. The above discussion suggests that C and hence \B V / B z \ at the disc surface is large. The actual 
magnitude of \B V \ is then controlled by the requirement of a force free equilibrium in the corona. This results in a significant 
departure of the poloidal field from the original stellar dipole through a significant inflation and opening out of the field lines. 

Equation (^) is solved for Rd < r < _R cxt where Rd and 7? cxt correspond to the inner and outer disc radii. Ra is usually 
taken to be small enough that the magnetic torque that is applied there overwelhms the viscous torque so that angular 
momentum transport is magnetically regulated. 



2.3 The location of the inner disc radius 

An accurate calculation of the inner disc radius would require the solution of the accretion disc structure problem with the 
inclusion of all magnetic stresses and allowing for sub-Keplerian rotation in the innermost disc region. In most studies an 
approximate value for Rd has been calculated based on the criterion of disc disruption at a radius where magnetic stresses start 
to dominate viscous stresses. In order for stellar accretion to proceed the magnetic field must not disrupt the disc exterior to 
the corotation radius since the magnetic torque for R> R c imparts angular momentum to the disc gas which cannot connect 
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to the stellar field. Assuming accretion takes place so that Rd < R c , the inner disc radius can be estimated by requiring that 
the magnetic torque balances the rate of angular momentum advection by the flow or that: 



M 



d (yiR 2 
Ir~ 



t[b v b z 



=i? d -Rd 



(16) 



where M is the mass accretion rate and the positive sign corresponds to clockwise rotation. Because of the rapid increase of 
the magnetic field strength inwards estimates of Rd using equation ( |l6| ) (e.g. Wang 1995| ; Yi 1995) give Rd close to R c when 
the star is near its equilibrium spin rate. 

However, we point out that this analysis is based on an assumed stellar dipole structure for B z . When the effect of an 
increasing B v at the disc surface is taken into account we expect that both B z and the magnetic torques will be modified. 
We will return to this point below. 

In our numerical calculations we have mostly set Rd = R c implicitly assuming that the star is close to its equ ilibr ium 
spin rate. For completeness we also performed some calculations with a sub-Keplerian inner disc region (see Section 4.3). 



3 NUMERICAL SOLUTION 



3.1 Previous work 



Equation (Ji]) has been solved by Bardou & Heyvaerts (1996) using a relaxation method with a coordinate transformation 
r ^ 1/r so that the solution can extend to arbitrary large radii in principle. Bardou & Heyvaerts (1996) adopted equation 
( p^| ) for Btp/Bz with a maximum expected value of C ~ 3.25 • 10 3 (Rd/R) 1 ^ 8 . However, in practice the numerical method 
would not converge unless C was multiplied by a small parameter A whose value depends on R c /Rd- As a result the maximum 
value of \B lp /B z \ at the disc surface that could be obtained in a converged calculation was 2.55. For this particular case 
A = 1.6 • 10~ 5 (A. Bardou, private communication). Bardou & Heyvaerts (199C) found that the force-free equilibria consist of 
closed field lines whose footpoints are shifted radially outwards w ith respect to their position in the dipolar state and with 
field lines that are flattened at small 9 as compared to dipolar ones. Bardou (1999) found the same indication from a similarity 
solution. Numerical convergence problems were encountered also by Wolfson fe Low (1992 ) who used a relaxation method 
to solve equation (^|) with a prescribed / in the context of the solar corona. They found that no equilibrium solution could 
be calculated for \B lp /B z \ <; 1. We comme nt that methods for solv ing equation (^|) based on a simple iteration method for 



producing a contraction mapping (see also Wolfson fe Verma 1991 ) using successive previous approximations for the right 
hand side tend to fail to converge if //' is a sufficiently rapidly varying function of 9. By comparison with solutions of related 
algebraic problems, the lack of convergence in solving equation ^ is not necessarily related to the emergence of additional 
solutions through a bifurcation. On the other hand studies that employ different methods of solution do not seem to encounter 
convergence problems (i.e. 



Klimchuck 1991 



floumeliotis, Sturrock & Antiochos 1994) 



We conclude that previous studies of force-free equilibria related to star-disc interactions were restricted to modest values 
of \B V /B Z \ at the disc surface due to restrictions enforced by the adopted numerical method. 



3.2 Numerical method 

In this paper we adopt a different method of solution of equation (^) to those indicated above. In our case equation (^|) can 
be solved without the use of a controlling parameter A and in principle for any disc surface value of \B V /B Z \. In practice some 
restrictions also apply in this case but we were able to perform calculations for discs with a ratio of inner to outer radius of 
several decades and for a maximum value of \B V /B Z \ ~ 10 2 . 

In our method equation (^) is transformed into a parabolic PDE by moving the source term to the left of equation ^ 
and by equating the new expression to a multiple of the time derivative of which is expected to tend to zero close to an 
equilibrium configuration. We use an explicit finite-difference scheme to integrate the parabolic PDE forward in time as an 
initial value problem until a steady state is achieved. Following this approach rather than solving the elliptic equation ^ 
directly, we were able to obtain equilibria corresponding to a wide range of values of the source term. 

The total magnetic stream function, hereafter denoted by can be written as a sum of contributions arising from disc 
currents, and the central dipole, ty d p- Thus *t = ^dp + *d, where ^d P is given by: 

*d P = -K^—^-, (17) 
r 

and A' is a normalizing constant. Note that A^d P = since / = for a dipole field. Equation ^ with / 7^ is then solved 
for *d only from: 

A* d = -/(*t)/'(*t). (18) 
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We find it convenient to drop the subscript "d" from the disc magnetic stream function and use dimensionless quantities, 
f — r/r*, = \&/\I/o an d B = B/Bq where B represents the magnitude of either the toroidal or the poloidal magnetic field 
components and we take = Bor^ such that Bo = Bd P (r*,0), the magnitude of the dipole field on the equatorial plane. 
Equations (jij) and (|l^) can then be read in dimensionless form if all quantities are replaced by their barred equivalents. For 
simplicity we will drop the bars from the dimensionless quantities from now on. 

The computational domain is bounded by an inner radius n n , an outer radius r out = -Rext, the disc midplane at /i = 
and the symmetry axis, fj, = 1. The grid spacing is taken to be uniform in both r and fi. 

At r = ri n and fi = 1 respectively we specify the boundary condition $ = 0. The first condition specifies only the dipole 
flux exists at the inner boundary while the second is a requirement of the coordinate system. On the disc midplane, /i = 0, 
we impose the symmetry condition, dty/dfi = 0. 

We have considered two different outer boundary conditions. The first one is the Dirichlet condition: *I> = which forces 
the field due to the disc to be tangential there. Physically this corresponds to a conducting boundary which excludes the field 
arising from currents in the disc but not the original dipole field which may be assumed to have had infinite time to diffuse. 
The second outer boundary condition that we use is the Neumman condition: 4-^ = 0. This makes the disc field radial at the 
outer boundary as might be the case if there were a coronal wind there. The second boundary condition leads to much more 
open configurations than the first and hence to larger departures from the stellar dipole field. 

We add a term Fftty/dr) to the right-hand side of the dimensionless form of equation ( |f8| ) where we are allowed to 
choose F as an arbitrary function of position and solve the equation with a forward in time, r, and centered in space (FTCS) 
finite-difference scheme on a (Nr x Nm), (r, p) computational grid. The finite-difference equation that we use is given at a 
grid point denoted with subscript by: 



'•». 



1,3 



dr (dx) 



where dr, and dfi, are the uniform grid spacings in r, and n, respectively. Variables without superscripts are at r level n. The 
time step is dr. We assume that the disc is truncated at Rd by an infinite conductor so we have B v = for 1 < R ^ Rd- For 
all other values of R on the equator (fj, = denoted by subscript 1) we use equation dli| ) approximated as follows: 

- B vli+l/2,l = C Bz \i+1/2,1 ' ( 20 ) 



In equation ( po| ) r^+1/2 is either given by the Keplerian rate or by the form proposed by |Campbell (1987 ) (see Section 4.3) 



fi* is calculated on the (i + 1/2, 1) grid point that is closest to the prescribed value of R c . The vertical magnetic field at the 
disc surface is calculated at the (i + 1/2, 1) points by numerical differentiation of equation (^). 

Using equation ( po) ) we calculate /i+1/2.1 = Ri+1/2 ^vli+1/2 1 while (ff')i,l is calculated from: 

1 /f+1/2,1 ~ /i-1/2,1 , > 

2 W *li+1/2,1 _ W t| l -1/2,1 

Since / is a function of "Jt only we construct a table of values of / as function of ^tL ! on the equator. We then use this table 
to calculate / for any required value of ^ t in equation ( p^j ) using linear interpolation. 

The main disadvantage of the method we use is the constraint on the time step coming from the requirements of numerical 
stability. We find that the maximum allowed timestep decreases rather steeply with C. 

In order to achieve a more rapid evolution for the same computational time and always keeping in mind that we are 
only interested in final steady state equilibria, we have introduced a spatially varying "diffusion coefficient", (1/F), that 
multiplies the right-hand side of equation (|l9|). This coefficient is taken to be equal to dTij / mia(dTij) and therefore allows 
for a more advanced evolution of the values at the grid points where dr^j, the maximum allowed timestep based on local 
stability considerations, can be larger. Since the source term peaks near the corotation radius but it is significantly smaller 
elsewhere it is possible to reach the equilibrium in significantly smaller computational time in this way. 

Time integrations have been performed until the right hand side of equation ( |li| ) becomes smaller than ~ 10 -6 in 
dimensionless units although the magnetic field configuration is usually found to have converged to its final form well before 
this criterion is satisfied. 



4 NUMERICAL RESULTS 

The parameters used in the calculations presented here are summarised in Table [j]. All the calculations were performed with 
ri n = 1. Different values were used for i? e xt and C and calculations were done using both outer boundary conditions on 9. 
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B.C. 


C 


Rcxt 




Nr, Nm 


B(p 1 B z | (i=o 


X st 


3 


11 


1 


200,80 


2.56 






21 






2.83 






60 












OU 






9 Q7 




10 


60 




100,50 


9.85 






80 






9.90 




30 


11 




50,30 


25.11 




30 


21 




120,50 


28.18 




30 


60 


0.875 


100 50 


29 60 




30 


11 


0.2 


50,30 


23.86 




30 






100,50 


24.90 




60 




1 


50,30 


50.21 




60 


60 


0.875 


100,50 


59.21 




120 


11 


1 


50,30 


100.42 




120 


60 


0.875 


120,50 


118.30 


2 nd 


3 


60 


1 


200,80 


2.96 




10 




0.875 


120,50 


9.86 




60 








59.15 




120 








118.30 



Table 1. Parameters used in the numerical calculations described in the text. All calculations above the third pair of horizontal lines 
were performed with the first (Dirichlet) outer boundary condition: <t = at r = Rext- Those below were performed with the second 
(Neumman) outer boundary condition which has the normal derivative of f on r = Rext vanishing. For both outer boundary conditions 
the maximum value of \B V /B Z \ in the calculations, given in the last column, is 118.3. Additional calculations to those listed here were 
performed in order to study the variation of the torque between the star and the disc with C. 



All the calculations with Q. — Qk were performed with Ra — R c . Our choice of R c — 3 implies tt* — 0.19S!k(-R = R*) which 
corresponds to a rotation period of ~ 5 days for R+ = 3Rq and Af* = 0.5M@. The observed rotation periods of CTTS are 
3-12 days (Bouvier, Forestini & Allain 1997 and references therein). 

The choice of 7? ex t was dictated mainly by the limitations of the computational method since a much larger value of 7? ex t 
would require increased values of Nr and Nm for the same resolution to be achieved as in the cases with smaller outer disc 
radii. 

We chose C = 3-120 in order to examine the effect of an extended range of values of \B lp /B z \ at the disc surface on the 
magnetospheric equilibrium. Various authors have assumed \B V /B Z \ ~ 1 at the disc surface in their a nalysis. |G!hosh fc Lamb 
(1979^) argue that anomalous resistivity in the disc would limit B v jB z . Aly & Kuijpers (199C) and Livio & Pringle (1992) 
assert that fast reconnection outside the disc would lead to a similar result. These physical arguments have not so far been 
supplemented by detailed calculations in the context of star-disc interactions. Our maximum value of C was chosen so that 
the computational times required are reasonable. Nevertheless much larger values of C and therefore of B v on the disc surface 
could result in very large toroidal magnetic pressure on the disc surface and the subsequent destabilisation of the disc. 

In Table [l] we have also included the parameters of calculations performed with a sub-Keplerian form of Q. This is 
characterised by the "fastness" parameter u = /fiK(-Rd) (for more details see Section 



4.1 Qualitative considerations of star— disc interactions 

The main result of the production of toroidal field in the star-disc corona is the presence of a feedback mechanism that reduces 
\B Z \ at the disc surface. In this way the magnitude of B v does not become too large to maintain a force free equilibrium when 
C is very large. As C is increased, the field becomes increasingly twisted. The result is the production of azimuthal currents 
as well as a toroidal field component in the corona. The azimuthal currents tend to reduce \B Z \ near the disc and hence the 
magnitude of B v . To see this in a qualitative way: the disc flux obeys equation ([H|) which, for the boundary conditions used, 
can be converted to an integral form: 

tf= / G(r, M ,r', Al ')/(* t )/'(* t )r' 2 d/i'dr'. (22) 
Jv 
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Figure 1. Results with C = 3. Left Panel: Radial profile of \B V \ at the disc surface with i? ex t = 11 (solid line), 21 (dashed line), 60 
(dot-dashed line) and 80 (dotted line). Right Panel: Radial profile of B z (r, 0) with i? cx t = 11 (solid line), 21 (dashed line), 60 (dot-dashed 
line) and 80 (triple dot-dashed line). The initial dipole profile of B z (r, 0) is represented by the dotted line. 

Here G is an appropriate Green's function and the stream function inside the integral, taken over the computational domain, 
is evaluated using the primed coordinates. Thus: 

*t = J G(r, fj,, r',fi')f(<f t )f'(^ t )r' 2 dfi'dr' + * dp . (23) 

For small C, f is small and is close to the dipole value 'I'dp- However, for large C and hence /, a solution can only be 
found (or the force free equilibrium can only be maintained) by otherwise reducing the magnitude of the integral. Note that 
if *d P is scaled by multiplying by a constant value, both and the integral above scale similarly. The magnitude of the 
integral can be reduced by having fewer field lines intersecting the disc and therefore more field lines for which / = 0. If the 
outer boundary condition does not allow that, the magnitude of the integral can also be reduced by the field lines intersecting 
the disc at progressively larger radii (and smaller field magnitudes), as C increases. In either case, the magnitude of B v is 
controlled and the poloidal field takes on either a more inflated or a partially open structure. 

We now discuss our results for a disc with Keplerian rotation everywhere and _Rd = i? c = 3 (those cases are annotated 
with u = 1 in Table [j]). 

4.2 Keplerian disc rotation profile 

The calculations presented here were performed using five values of C in the range 3-120 and using the Dirichlet outer 
boundary condition. The maximum value of J? cxt used in these calculations decreases with C (see Table [i]). As C increases 
the resolution required to resolve fully the magnetic field structure around the corotation radius increases. We have therefore 
restricted the external radius for C = 30 or larger so that a reasonable resolution of the area around corotation can be 
achieved. 

The maximum value of \B V /B Z \ at the disc surface found in the calculations with a Keplerian disc rotation profile is 
~ 100 being much larger than values considered by previous authors. Consequently we study the effect of increasing \B Lp /B z \ 
on the force-free magnetospheric field for a large region of parameter space. 

Below we present results obtained with C — 3 corresponding to maximum value of \B V /B Z \ — 2.56-2.97 for _R oxt = 11-80. 
When C = 3 we were able to achieve the highest resolution and the largest range of i? cx t. 

4-2.1 Results with C = 3. 

All the calculations with C = 3 were performed with Nr = 200 and Nm = 80 therefore cases with smaller i? oxt have better 
resolution. 
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Figure 2. Isocontours of in the R = r(l — /i 2 ) 1 / 2 , 2: = r^i plane. The contour levels correspond to initial values of >Pt (dipole contours 
plotted with dashed lines) for specific values of r (and fi = 0). Those are: Upper left panel: r = 1,2,4,6,8,10,11, Upper right panel: 
r = 1, 4, 6, 8, 10, 15, 21, Lower left panel: r = 1, 4, 6, 8, 15, 40, 60, Lower right panel: r = 1, 4, 6, 8, 15, 40, 80. 



In the left panel of Fig. |l| we plot the equilibrium radial profiles of \B V \ at the disc surface. Since these profiles are well 
resolved we can conclude that increasing the computational domain leads to smaller values of \B V \. This follows from the 
dependence of B z on 7? ox t . In the right hand panel of Fig. [l] we plot the radial profile of B z at the disc surface for the different 
values of R cx t- 

The poloidal field is similar to the stellar dipole field inside corotation where B v — but it deviates from that significantly 
at all other radii. Immediately outside corotation the field becomes smaller than the stellar dipole field as a result of field line 
inflation. At the largest radii the poloidal field eventually becomes larger than the dipole field locally. The radius at which 
this transition occurs moves to larger radii as i? C xt increases indicating the influence of the outer boundary condition. The 
structure of the magnetic field is not radially self-similar. The departure from self-similarity in the inner part of the disc is 
due to the transition from B v = inside the corotation radius to B v given by equation ([HJ) for R > R c while the outer part 
of the disc is influenced by the outer boundary condition. 

Field line expansion as seen in Fig. (H) for all R ext used is a consequence of field line twisting together with the assumption 
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C = 3 




R. 

Figure 3. Rf/Ri as a function of Ri, where Rf is the force-free and Ri the dipole footpoint radii for C = 3 with i? cx t = 11 (open circles), 
21 (filled circles), 60 (open triangles) and 80 (filled triangles). 



of a force-free equilibrium (see our discussion in Section 4.1). The inflation of field lines increases with increasing C or field 
line footpoint shear until the field lines become partially or fully open (Aly 1985; Wolfson & Low 1992; Newman, Newman 
& Lovelace 1992; Lynden-Bell & Boily 1994; Aly 1995; Wolfson 1995). In the fully open state line twisting disappears and 
B v = 0. 

In Fig. ^| we plot the ratio of the final Rf to the initial or dipole Ri footpoint radius as a function of Ri for a group of 
field lines. The values of Ri chosen are the same as in Fig. ^ with Ri ^ 4. Different curves correspond to different values of 
Tiext. The ratio Rf/Ri increases linearly with Ri until a turnover radius is reached where the outer boundary starts affecting 
the rate of inflation. For Ri larger than the turnover radius Rf/Ri decreases with Ri. For the cases with i? C xt = 60 and 80 the 
values of Rf / Ri are similar for Ri ^ 10. We therefore conclude that the maximum inflation has been reached for the innermost 
disc region in these cases. The maximum overall value of Rf/Ri is 3.2 and it corresponds to \B tp /B z \ of 2.73. The maximum 
value of Rf/Ri quoted by Bardou & Heyvaerts (1996) is 1.9 but this is an underestimate because of resolution problems 
(A.Bardou, private communication). Due to the effect of the outer boundary in our calculations it is difficult to predict how 
inflated the configuration would become for R ext — > oo. However, we expect the field lines will remain closed within a radius 
of approximately 3R C for C — 3. 



4-2.2 Results for larger values of C 

We performed calculations for C = 10 with i? B xt = 60 and 80. Both these have Nr — 100 and Nm = 50 so the resolution is 
reduced with respect to the C — 3 case. In Fig. ^ we plot the radial profiles of ty t {r,0) and B z (r, 0) for C = 10 and also for 
C = 3 and R cxt = 80. For i? cx t = 80 we find that 9t(r, 0) is smaller for C = 10 than for C — 3 everywhere in the disc as 
expected. For the smaller 7? ex t value the effect of the outer boundary is stronger as it was found also for C = 3. We find that 
B z decreases with C for R < 15. For larger radii the trend is reversed because of the effect of the outer boundary. 

In the left panel of Fig. ^ we plot dipole and equilibrium poloidal field lines for C = 10 with _R cxt = 80. Field lines interior 
to 2R ext /3 are characterised by the largest flattening at small 8 while the outer field lines are more spherical due to the 
effects of the outer boundary. A value of i? cx t exceeding the outer disc radius by a factor 10-100 might be necessary in order 



to remove the effects of the outer boundary on the calculation of field line inflation (see Roumeliotis, Sturrock & Antiochos 



1994) 



In the right panel of Fig. |5j we plot Rf/Ri as function of Ri for C = 3 and C = 10, each with i? cx t = 60 and 80. For Ri ^ 8 
the ratio Rf/Ri is 1.75 times larger for C = 10 than for C = 3 showing greater inflation. But for larger values of Ri this factor 
drops to 1.25. The larger influence of the outer boundary on the equilibrium for C = 10 is manifested in the smaller turnover 
radius observed in this case. It is not clear if the ratio of final to initial footpoint radius found for C = 10 at _Ri = 6 would 



continue to increase linearly with Ri if the outer boundary effect was not present. In the Bardou & Heyvaerts (1996) results 
Rf/Ri increases linearly with Ri for R^ R c . In our case \B V /B Z \ tends to a constant for R 3> Rc and Rf/Ri could behave 
similarly. 
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Figure 4. Left Panel: Radial profile of ^t(r, 0) for C = 10 with -R cx t = 60 (solid line) and 80 (dashed line) and for C = 3 with B ex t = 80 
(dot-dashed line). Right Panel: Radial profile of B z (r, 0) for the same parameters as in the left panel. In both panels the initial dipole 
profiles are represented by the dotted line. 




Figure 5. Left panel: Isocontours of \I/t on the R, z Diane plotted for C = 10 with i? ex t = 80. The initial dipole contours are plotted 
with dashed lines. The radii chosen are given in Fig. El Right Panel: Rf/Ri as a function of R\ for C = 3 with i? ex t = 60 (open circles) 
and 80 (filled circles) and for C = 10 with i? ex t = 60 (open triangles) and 80 (filled triangles). 



We now discuss results obtained for C = 30-120. These calculations were performed with Nr = 50 and Nr = 30 and 
Rext = 11 except for one case with C = 30 and J? ex t = 21 which had Nr = 120 and Nm — 50. The disc surface equilibrium 
profiles of B v for C = 3-120 and J? oxt = 11 are compared in Fig. |(| Note that the magnitude of B v varies only by a factor of 
3 as C varies between 3 and 120. 

The effect of the growing \B V \ on the stream function for C ^ 30 and for i? cxt = 11 is to make it larger than the 
dipole value, more so around R = R c , leading to an enhancement of B z (r, 0) there. The disc surface radial profiles of B z for 
C — 3-120 are compared in the left panel of Fig. 0. Because the toroidal magnetic pressure increases sharply just outside 
R — R c it is expected that a larger poloidal field just interior to R c is required for equilibrium. Thus field line deflation results. 
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Figure 6. B v as a function of radius for C = 3 (solid line), C = 30 (dashed line), C = 60 (dot-dashed line) and C = 120 (triple 
dot-dashed line). 
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Figure 7. Radial profile of B z (r, 0). Left panel: C = 30 (solid line), C = 60 (dashed line), C = 120 (dot-dashed line) and C = 3 (triple 
dot-dashed line). All cases have -Roxt = 11. Right panel: C = 30 with R e xt = 11 (solid line) and with R cx t = 21 (dashed line) and C = 3 
(dot-dashed line). The initial dipole profile is represented by the dotted line. 



This occurs in this model because the field lines are not sheared inside corotation. For R > R c B z (r, 0) becomes smaller than 
the dipole value before increasing close to the outer boundary. Due to the increased \B V /B Z \ at the disc surface, B z decreases 
with C outside corotation as expected. 

To study the dependence on the outer boundary radius for large C, we performed a calculation with C — 30 and i? cx t — 21 
and compared the result to the case where R ex t = 11. The radial resolution is approximately the same for these cases (see 
Table [I]) . The corresponding radial profiles of B z (r, 0) are plotted in the right panel of figure 0. When _R C xt is increased B z (r, 0) 
values are unmodified inside R = 10 so the deflation of the field lines there is independent of the outer disc radius. The exterior 
point at which B z (r, 0) exceeds the dipole value moves to a larger radius for C = 30 as compared to that when C = 3. There 
is also a change to the radial dependence of B z (r, 0) which becomes flatter at large R for C ^ 30. 
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In Fig. ^]we plot the magnetic field lines for C = 30 with R cxt = 11 (left panel) and i? cxt = 21 (right panel). We see that 
the deflation of field lines for r < 10 is the same in these cases. At larger radii inflation rather than deflation is observed for 
J?ext = 21. The field line inflation observed at r > 10 for C — 30 with i? cxt = 21 is smaller than that found for C = 3 and 
the same i? cxt (compare with Fig. ^). This result is unexpected given the larger \B ip /B z \ at the disc surface in the C = 30 
case. This is probably the result of the increased influence of the outer boundary condition on the equilibrium configuration 
for increasing C values. 

In the next section we will discuss results obtained using a sub-Keplerian form of f2 and which results to a smoother 
radial variation of B v / B z on the disc surface. 



4.3 Sub— Keplerian disc rotation profile 

As the magnetic field strength increases with decreasing disc radius and the magnetic stresses start to affect the radial and 
vertical disc equilibrium the disc is expected to ultimately corotate with the star. Interior to the corotation radius the disc will 
then be sub-Keplerian. In the present study we investigate the effect of sub-Keplerian rotation on B^/Bz and the equilibrium 
poloidal magnetic field. We modify fi so that it becomes sub-Keplerian for R ~ R c . The form of fi that we will use is that of 



Campbell (1987) which is given, for R > Rd, by: 

. / R Nil 

(24) 



n = n K (R d ) | (^) 3/2 - (i - ^)ex P - c)- 1 (JL - 1)] J 



Here lu = Q*/Q,K(Rd). At R = Rd we have Q = and dQ,/dR — 0. Interior to R — Rd, the material is taken to rotate with 
fi = f2*. We plot Q, as a function of R in the left panel of Fig. ^ for R c = 3 and to — 0.875 and 0.2. The condition: Q, = 0,+ at 
R — R c determines Rd through the chosen value of w. 

For the largest value of ui, Rd almost coincides with R c and fi remains Keplerian until very close to R c . For smaller values 
of w we have Rd < R c and in the case where lo = 0.2, Rd almost reaches the stellar surface (R — 1). Here as in the previous 
sections we have B v — for R < Rd but a part of the differentially rotating disc may exist inside corotation and that can 
affect the equilibrium magnetic field. The variation of B v /B z with r at the disc surface can be seen in the right panel of 
Fig. where we took C = 30. For u = 0.875 the profile of B v /B z is virtually unchanged with respect to the Keplerian case 
with Rd = -Rc- A larger effect on B v jB z occurs when u> = 0.2. In that case there is a region inside corotation where B v /B z 
reverses sign. 



4.3.1 Results with uj = 0.875 and C = 30-120 

These calculations have Rext = 60. In the left panel of Fig. ^| we plot B z (r, 0) for C = 
© 0000 RAS, MNRAS 000, 000-000 
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Figure 9. Radial profiles of Q (left panel) given by equation (M) and of corresponding B v / B z (right panel) using R c = 3 with ui = 0.875 
(dashed line) and w = 0.2 (dot-dashed line). The radial profiles of S!k and corresponding B v /B z are plotted with a solid line. 
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Figure 10. Radial profile of B z (r, 0) (left panel) and of <f t (r, 0) (right panel) for u> = 0.875 and C = 120 (solid line), C = 60 (dashed 
line), C = 30 (dot-dashed line) and C = 3 for Q = !!k (triple dot-dashed line). The initial dipolc profiles are represented by a dotted 
line. 



with the corresponding ^t(r, 0) profiles shown in the right panel. We also plot the case C = 3 with Q = Qk- Apart from 
boundary effects, \B Z \ decreases with C, the decrease being faster between C = 60 and 120 than between C = 30 and 60. 
These results reinforce our expectation: B z — > for C — > oo. 

As seen in the right panel of Fig. [H^ the disc values of $ t increase with C between C — 30 and 120. This is contrary 
to theoretical prediction and to the trend observed between C — 3 and 30. As we will see in Section [T4| this result depends 
strongly on the outer boundary condition. 

The profiles obtained in the non-Keplerian case with u = 0.875 are smoother around corotation than in the Keplerian 
case. This is due to the smoother radial variation of B v /B z . However, B z (r, 0) is still larger than the dipole value there and 
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Figure 11. Left Panel: Radial profile of B v at the disc surface for C = 30 with B ex t = 11 and for u> = 0.2 with JVr = 100 (solid line) 
and with Nr = 50 (dashed line) and for f2 = Ok (dot-dashed line). Right Panel: Radial profile of B z (r, 0) for the same parameters as in 
the left panel. The initial dipole profile is represented by the dotted line. 

some deflation is observed interior to r — R c . At larger radii field lines are inflated as we would expect given the increase in 

Rcxt- 

4.3.2 w = 0.2 

In this case Ra/R c = 0.34 and an inversion of the sign of B v may significantly affect the variation of B z with radius. We 
adopted C = 3 with R cxt = 21 and C = 30 with J? cxt = 11 with two different resolutions (see Table [l] for details). For C = 3 
the difference between Keplerian and non-Keplerian cases is not significant. 

In the left panel of Fig. [n] we plot the disc surface profiles of B v for C = 30 and 10 = 0.2 obtained with two different 
resolutions and for C = 30 with fl = S1k- For C = 30 the larger radial variation of \B V \ in the non-Keplerian case has a more 
significant effect on the equilibrium. As the radial gradient of B^ increases for small radii a small deflation of the field lines 
is observed there. Eventually B^, decreases with radius and inflation is observed for R > 2. The reverse was observed in the 
Keplerian case where field lines were found to be deflated everywhere in the disc (always for 7? ex t = 11). 

In the right panel of Fig. [H] we plot the radial profile of B z (r, 0) for the same parameters as in the left panel. We note 
that when u = 0.2 the field deviates more from the dipole form around corotation, becoming larger than the dipole value 
locally, than in the case with f2 = Qk- This difference is mainly due to the larger radial variation of \B V \ with radius in the 
non-Keplerian case. The difference observed between the profiles obtained in the high and low resolution cases is mostly due 
to the radial shift of the gridpoint corresponding to R = Rc. 

We conclude that the deflation observed in previous calculations with a Keplerian rotation profile and with u> = 0.875 is 
related to the restricted twisting of field lines inside the corotation radius. In the case with u = 0.2 we have B v 7^ almost 
everywhere in the disc. Consequently field lines are much more inflated in this case. 

4.4 A different outer boundary condition 

To study the effect of the outer boundary condition on our results we have also performed calculations with the Neumann 
condition: = and for C = 3, 10, 60 and 120. The parameters used in these calculations are given in Table ^ All cases 
have ui = 0.875 except when C = 3 where a Keplerian rotation profile was used. The new boundary condition makes the disc 
field normal to the outer boundary. Thus we expect some of the field lines to be open at the equilibrium, at least at large disc 
radii far from the disc midplane. 

Indeed the results show much larger inflation and opening of the field lines in this case than in the cases presented in 
previous sections where the Dirichlet condition was used. Also almost no deflation is observed inside corotation showing the 
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Figure 12. *t(r, 0) (left panel) and B z (r,0) (right panel) for C = 3 (solid line), 10 (dashed line), 60 (dot-dashed line) and 120 (triple 
dot-dashed line). Note that the ^t{r, 0) profiles for the last two cases approximately coincide. The initial dipole profiles are represented 
by the dotted line. 



global effect of the outer boundary condition. As seen in Fig. |l2] both $t (r, 0) and B z (r, 0) decrease with C monotonically. 
For C = 120 the equilibrium seems to have reached a limiting configuration with almost all field lines open. 

Comparing the disc radial profiles of *I/t (r, 0) depicted here with those in Fig. [j^ we note that $ t {r, 0) decreases consider- 
ably when the second boundary condition is used. Also $t(r, 0) becomes approximately independent of r everywhere outside 
the corotation area which is the limiting behaviour that we expect for large values of C. Comparing the disc radial profiles 
of B z (r, 0) between figures and we note that for C = 3 the magnetic field is similar for the two boundary conditions 
for r ^ 10. Further out B z (r, 0) is smaller for the second boundary condition. This difference increases with C. We also note 
that B z (r, 0) has a stronger dependence on r when the second boundary condition is used. 

In the first three panels of Fig. |l^ we plot the magnetic field lines for C = 3, 10 and 60 respectively, using the same initial 
(dipole) field lines as in Fig. ^ (with R cx t = 60). For C = 3 the level of inflation is comparable for the two boundary conditions 
for r 10. For larger initial footpoint radii field lines become much more inflated for the second outer boundary condition and 
eventually they open up for r ^ 15. As C increases more field lines become open as seen by comparing the cases with C — 10 
and 60. The case with C = 120 is virtually identical to the C — 60 case. Some field lines remain closed even when C = 60 
as can be seen by plotting more contour levels (see lower right panel of Fig. |l3|). The closed field lines correspond to initial 
radii which are very close to the inner disc radius. Consequently //' is non-zero for an increasingly smaller range of values 
as C increases. Calculations for larger values of C could only be performed accurately using considerably larger numerical 
resolution. However, our present calculations have already shown that field lines at equilibrium will be approximately open 
for \B V /Bz \ ~ 10 2 if a Neumman outer boundary condition is used. 

The ratio \B V \/B P is found to be less than unity for /i < 0.1. Therefore B v m in most of the corona as expected on open 
field lines. In the limit of a large twist the disc field acts as to expel the poloidal magnetic field from the disc and consequently 
to annihilate the toroidal magnetic field, as long as the disc field is permitted to penetrate the outer boundary. In the limit 
of large C the poloidal field resembles the solution Z\ given by Low (198£). This corresponds to a perfectly conducting disc 
with central hole immersed in a central dipole field where some flux escapes to infinity such that the field is non singular at 
the disc inner edge. 



4.5 Modification of the spin— down and spin— up torques 

For the force-free equilibria calculated above, the poloidal magnetic field differs significantly from the original stellar dipole 
form. As a consequence, the local magnetic torque acting on the disc may also differ significantly from what one obtains 
assuming B z ~ Bd P - Since this torque regulates the spin evolution of the central star it is important to examine the possible 
effect of our results on the calculation of the total torque experienced by the star. 

The total magnetic torque JV can be split into two components of opposite sign. The first arises inside corotation where 
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Figure 13. Upper left to lower left panels: Isocontours of \Pt on the R, z plane plotted for C = 3, 10 and 60 using the initial footpoint 
radii given in Fig. H. Lower right panel: Isocontours of \I/t for C = 60 using 20 equally spaced contour levels. 



angular momentum is transferee! from the disc to spin up the star. The second arises outside corotation where angular 
momentum is transferred to the disc. In order for stellar accretion to occur we require < R c . In one set of calculations 
we assumed = R c and adopted Keplerian rotation. In that case the spin-up torque would be zero. In reality we expect 
Rd/Rc ~ 0.91-0.97 (i.e. Wang 1995) giving rise to a relatively small spin-up torque. This torque is probably small compared 
to that arising from direct accretion onto the star. 

We now calculate the total torques for the case where the disc is assumed to be Keplerian and truncated at the corotation 
radius. 



4- 5.1 The Keplerian case 



The total magnetic torque on the star is in general given by (e.g. Shosh & Lamb 1979b): 
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N = - B lp B z R 2 dR. (25) 

We adopt dimensionless parameters such that the torque is given in units of N* =BM- Of course all field values are calculated 
on the disc surface (fj, = 0). For the cases presented here the first boundary condition has been used. When = R c , as is 
considered here, N is positive definite and only acts to spin down the star. 

The total torque N can be calculated analytically for a dipole field with B v given by equation (jlij) and is given by: 

- ' ' 1 ■ R~\ (26) 



jV dp = — 
3 



2 / itcxt \ / Jlcxt \ 

3V~rT/ ~ V~rT/ 



In the left hand panel of Fig. |T| we plot N/C and N dp /C as a function of C for i? cxt = 11. Although N dp /C is constant 
with C, N/C decreases with increasing C as the poloidal magnetic field in the force-free equilibria decreases with increasing 
C. The ratio of force-free to dipole torque N/N dp for C ^ 30 varies between 0.14-0.013 which is significantly smaller than 
the same ratio in the case with C = 3. Although the effect of the outer boundary condition is to suppress inflation in all cases 
with large C values and _R ext = 11, the effect of the twisting of the field lines on the final equilibrium is an increased reduction 
of the total torque with respect to its corresponding dipolar value. 

In the right panel of Fig. [l4] we plot N and iV dp where we see that the behaviour of N as function of C is very different 
(in fact almost reversed) to that expected for a dipolar poloidal field. From this figure we see that N increases slowly for 
C 15 and it decreases for larger values of C. As the values of N for C ^ 30 are between 10 and 100 times smaller than the 
dipolar values the increasing field line twisting has an important effect on the spin down of the star. The magnetic breaking 
timescales calculated are longer and therefore magnetic breaking models might have to be reevaluated. However, we have not 
taken into account here the disc region inside corotation. Also the results described here for C ^ 30 have 7? ox t = 11 which is 
relatively small and therefore the effect of the outer boundary condition is important. 



4-5.2 The sub-Keplerian case 

In the case where Q ^ Qk the total magnetic torque includes both spin-up and spin-down components. For the case with 
lu — 0.875 the spin-up component is of no significance. In the following we will compare the total magnetic torques calculated 
using a sub-Keplerian rotation profile with lo = 0.875 to the corresponding dipolar values. Calculations with similar parameters 
were performed with both outer boundary conditions. 

In Fig. |l| we plot N/C and N dp /C as a function of C. The values in the left panel of this figure are from calculations 
performed with the first outer boundary condition. Here, as for the Keplerian case, we find that N/C decreases with C and 
that the ratio N/N dp ranges from 0.56 for C — 3 to 0.05 for C = 120. The larger ratios here are mainly a consequence of the 
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Figure 15. N/C (solid line) and N dp /C (dashed line) as functions of C. All cases have ui = 0.875 and i? e xt = 60. Left Panel: Results 
with the first outer boundary condition. The case with C = 120 has Nr = 120 while the rest of the cases have Nr = 100. Right Panel: 
Results with the second outer boundary condition. Many more models are employed here than those tabulated in Table [l]. The torque 
is seen to vary smoothly with C. 



use of a sub-Keplerian rotation profile. However, For C ^ 30 we have again N/C oc 1/C approximately. Therefore, the trend 
for a smaller total torque is clear and the difference between force-free and dipole cases is again significant for C> 1. For 
B<f, ~ B z the difference is not as large although we have to reserve final judgement until calculations with much larger values 
of i?ext are performed. 

For calculations performed with the second outer boundary condition we found significantly smaller values of B z (r, 0) 
which decreases faster with r than when the first boundary condition is used. The steeper decline of B z (r, 0) with r is more 
noticeable in the cases with large C values. For large C values we therefore expect much smaller values of N/C to arise when 
the second boundary condition i s us ed. In the right panel of Fig. |l| we plot N/C and N dp /C as a function of C for the 
calculations presented in Section 4.4. The ratio N/N dp varies in the range 0.45-0.003. Although this ratio is almost unity for 
C — 3, as was the case for the first boundary condition, it becomes progressively smaller for larger C values. For large values 
of C we have N/C ~ 1/C 3//2 when the second boundary condition is used. 

In summary, when \B V /B Z \ on the disc surface is large the magnetic torque acting on it will be much smaller than that 
estimated assuming B z ~ Bd P , especially when external conditions enable the opening of field lines. 



5 DISCUSSION 

We first summarise our results. For a Keplerian rotation profile with the first boundary condition the maximum inflation 



observed in our calculations with C = 3 is such that Rf/Ri = 3.2. This is a confirmation of the results of Bardou & Heyvaert 



(1996). The expansion of the field lines is restricted by the outer boundary in this case. However, for \B V /B Z \ ~ 3 at the disc 
surface the field lines with r ^ 37? c do not vary much once R cx t > 10. Thus although significantly inflated they will remain 
closed in the force-free equilibrium. Inflation increases with C so that when \B lp /B z \ — 10, a maximum Rf/Ri ~ 5 is found 
for Ri = 8. 

In calculations with C = 30-120 and 7? ex t = 11 a new effect was found. Inner field lines become deflated rather than 
inflated with respect to the initial dipolar configuration. We expect the level of deflation, and consequent increase in poloidal 
magnetic pressure, to depend on the rate of increase of toroidal magnetic pressure, which it has to balance around R — R c . 
We demonstrated, by varying 7? ext , that the deflation is not related to the outer boundary. We found that the level of field line 
deflation for r < 10 was restricted when a smoother non Keplerian rotation profile was used. However, some level of deflation 
could be typical of the field structure in the region near the corotation radius. That applies especially to discs with r]/v ~ 1 
and with Ra ~ R c when a Dirichlet outer boundary condition is used. 

In all the calculations presented here the equilibrium poloidal field is more nearly radial above the disc surface than the 
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stellar dipole field. The lengthening of the field lines is accompanied by a reduction of \B V \/B P for increasing [i. We expect 
that the force-free condition constrains \B V \/B V to remain of order unity even for \B tp /B z \ S> 1 at the disc surface. 

We also performed calculations with a Neumman outer boundary condition that required the disc field to be radial at 
the outer boundary as it might occur if a wind existed. In these cases, other parameters being equal, the equilibrium field 
structure was more open. The ratio \B v \/B-p was found to be less than unity for /i < 0.1 and w in most of the corona as 
expected on open field lines. In the limit of a large twist the poloidal magnetic field tends to be expelled from the disc and 
consequently the toroidal magnetic field is annihilated, as long as the disc field is permitted to penetrate the outer boundary. 

We also found that for large values of \B V /B Z \ on the disc surface the total magnetic torque is much smaller than that 
estimated assuming B z ~ Bd P - For C = 120 the torque is only ~ 0.003 of the dipole value when the second boundary condition 
is used and ~ 0.013 — 0.05 times the dipole value when the first boundary condition is used (depending on the form of fi). 
Since reasonable values of r\ for circumstellar discs tend to produce very large values of \B V /B Z \ on the disc surface, we believe 
that the modification of the total torque discussed here will be significant if a force free equilibrium can be attained. We 
therefore expect significant modification of the magnetic breaking times of T Tauri stars to result from the large twisting of 
the magnetospheric field lines. 



5.1 The spin— down of neutron stars 



Recent observations of accreting neutron stars ( Nelson et al. 1997 ; Chakrabatry et al. 1997 ) suggest that the total torque 
between star and disc can osc illate producing alternating s pin-up and spin-down phases with little evidence of any correlation 
with the mass accretion rate. Li fc Wickramasinghe (199c| ) have suggested that this may be due to variations in the structure 



of the magnetosphere. 

The work presented here indicates that a large range in the value of the spin-down torque may be obtained depending on 
the outer boundary condition and effective disc resistivity. In cases favouring large field line inflation and open field lines the 
spin-down torque resulting from the disc-magnetosphere interaction becomes very small. Thus if the magnetosphere oscillates 
between such a state and one with smaller inflation, oscillations in the spin-down torque and hence the total torque may be 
produced. The change from open to closed field lines may also affect the magnetic torque that arises from a wind, if such a 
wind is produced. Changes in the magnetospheric configuration might occur through reconnection through the disc midplane, 
variations in outer conditions being more or less favourable for open field configurations, or variations in the internal disc 
resistivity. 



ACKNOWLEDGMENTS 

This work was supported by the European Union grant ERBFMRX-CT98-0195. V.A. acknowledges support by the State 
Scholarships Foundation (IKY) of the Hellenic Republic through a postgraduate studentship. The authors are grateful to 
Anne Bardou for useful discussions. 



REFERENCES 

Aly J. J., 1984, ApJ, 283, 349 

Aly J. J., 1985, A&A, 143, 19 

Aly J. J., 1991, ApJ, 375, 61 

Aly J. J., 1995, ApJ, 439, 63 

Aly J. J., Kuijpers J., 1990, A&A, 227, 473 

Armitage P. J., Clarke C. J., 1996, MNRAS, 280, 458 

Bardou A., 1999, MNRAS, 306, 669 

Bardou A., Heyvaerts J., 1996, A&A, 307, 1009 

Bouvier J., Forestini M., Allain S., 1997, A&A, 326, 1023 

Cameron A. C, Campbell C. C, 1993, MNRAS, 274, 309 

Campbell C. C, 1987, MNRAS, 229, 405 

Chakrabatry D., Bildsten L., Grunsfeld J. M., Koh D. T., Prince T. A., Vaughan B. A., 1997, ApJ, 474, 414 

Ghosh P., 1995, MNRAS, 272, 763 

Ghosh P., Lamb F. K., 1979a, ApJ, 232, 259 

Ghosh P., Lamb F. K., 1979b, ApJ, 234, 296 

Klimchuck J. A., 1991, ApJ, 354, 745 

Li J., Wickramasinghe D. T., 1998, MNRAS, 300, 1015 

Livio M., Pringle J., 1992, MNRAS, 275, 244 

Low B. C, 1986, ApJ, 310, 953 

Lyndon-Bell D., Boily C, 1994, MNRAS, 267, 146 

Mikic Z., Linker J. A., 1994, ApJ, 430, 898 



© 0000 RAS, MNRAS 000, 000-000 



Accretion disc-stellar magneto sphere interaction 



Nelson R. W. et al., 1997, ApJ, 488, 117 

Newman W. I., Newman A. L., Lovelace R. V. E., 1992, ApJ, 392, 622 

Roumcliotis G., Sturrock P. A., Antiochos S. K., 1994, ApJ, 423, 847 

Sturrock P. A., 1991, ApJ, 380, 655 

Wang Y.-M., 1987, A&A, 183, 257 

Wang Y.-M., 1995, ApJ, 449, 153 

Wolfson R., 1995, ApJ, 443, 810 

Wolfson R., Low B. C, 1992, ApJ, 391, 353 

Wolfson R., Verma R., 1991, ApJ, 375, 254 

Yi I., 1994, ApJ, 428, 760 

Yi I., 1995, ApJ, 442, 768 



© 0000 RAS, MNRAS 000, 000-000 



